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Abstract.  Progress  in  computer  hardware  and  improvement  of  numerical  methods  made  solution  of  the  Boltzmann 
equation  for  rather  complex  gas  dynamic  problems  real.  The  method  developed  by  the  author  is  based  on  a  projection 
technique  for  evaluation  of  the  collision  operator.  The  computed  collision  integral  is  conservative  by  density,  impulse, 
and  energy,  and  became  equal  to  zero  when  the  solution  has  a  form  of  the  Maxwellian  distribution.  The  later  feature 
sharply  increases  its  efficiency,  especially  for  the  near  equilibrium  flows.  The  method  is  extended  on  a  mixture  of  gases 
and  the  gases  with  internal  degrees  of  freedom,  where  it  can  incorporate  real  physical  parameters  of  molecular  potential 
and  of  internal  energy  spectrum.  Examples  of  computations  for  a  range  of  Mach  and  Knudsen  numbers  are  presented. 


INTRODUCTION 

Solution  of  the  Boltzmann  equation  was  initiated  by  Nordsieck’s  group  in  Illinois,  USA  nearly  40  years  ago 
when  first  big  computers  have  appeared  [1].  In  this  method  velocity  space  was  divided  by  equal  cells,  and  collision 
integrals  were  evaluated  by  simulating  of  molecular  collisions  with  randomly  chosen  velocity  vectors.  Then 
computed  mean  values  per  a  cell  were  attributed  to  the  middle  points  of  the  cells,  and  the  obtained  system  of 
equations  was  solved  by  an  iterative  method.  A  few  years  later  another  method  in  which  the  collision  integrals  were 
evaluated  directly  at  the  nodes  of  the  grid  in  the  velocity  space  by  Monte  Carlo  technique  was  developed  by  the 
author  of  this  paper  [2].  Very  soon  the  main  deficiency  of  both  approaches  -  non-fulfillment  of  conservation  laws  in 
the  computed  integrals  -  became  evident.  In  [3],  a  splitting  finite-difference  scheme  for  the  kinetic  equation  has  been 
proposed,  and  then  a  special  correction  was  developed  to  satisfy  the  conservation  laws  at  the  relaxation  stage  [4], 
The  new  method  showed  itself  much  more  efficient  then  that  of  [2],  but  the  correction  introduced  some  additional 
numerical  viscosity,  and  required  artificial  assumptions  in  the  case  of  gas  mixtures  [5,  6].  The  next  step  in  the 
development  of  conservative  methods  has  been  made  after  construction  of  a  Discrete  Velocity  Model  [7].  In  this 
model  molecular  velocities  are  given  at  a  uniform  Cartesian  lattice  and  impact  parameters  of  a  molecular  collision 
are  such  that  post  collision  velocities  belong  to  the  same  grid.  Therefore,  it  imposes  a  restriction  on  impact 
parameters,  which  are  independent  in  the  true  Boltzmann  equation.  Based  on  a  similar  idea,  conservative  methods 
for  the  Boltzmann  collision  integrals  have  been  proposed  [8,  9]. 

The  conservative  method  without  any  selection  of  impact  parameters  has  been  developed  by  the  author  [10,  11]. 
It  was  then  improved  and  adopted  for  computing  near  equilibrium  flows  [12,  13],  extended  to  gas  mixtures  [14],  and 
gases  with  internal  degrees  of  freedom  [15].  The  main  features  of  the  method  are  described  bellow. 


NUMERICAL  METHOD 


Consider  the  Boltzmann  equation 

df  / dt  +  % -df  / dx  -  1 

It  is  solved  at  a  grid  of  N0  equidistant  nodes  £  with  a  step  h  confined  in  a  domain  Q  of  a  volume  V  .  On  the 
basis  of  Dirac’s  5  functions,  the  distribution  function  and  the  collision  integral  can  be  presented  in  a  form 


/(£> X, 0  =  X  fy ( x >  0<? (£  -4r)>  *>0  =  2  Jr 0£(£  -  ) 

y=l  y=\ 
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In  configuration  space,  an  arbitrary  discrete  grid  Xi  can  be  applied.  When  the  coefficients  of  expansions  (1)  are 
determined,  the  problem  reduces  to  the  solution  by  a  finite-difference  method  of  a  set  of  equations 

dfr/8t  +  ^r-dfr/dx  =  Ir  (2) 

Evaluation  Of  The  Collision  Integral 

Let  us  write  the  Boltzmann  collision  integral  for  a  mono  atomic  gas,  omitting  the  variables  X  and  t ,  in  a  form 

co  2/r  bm 

m  =  im®\=\  J  j  (f'f-'-ff-)gbdbd<pd& 

-oo  0  0 

where  velocities  after  collision  E, d '  are  related  to  those  before  the  collision  by  the  formulae 


d  =  €  +  k(k-g),d'  =  d~k(k-g),  g  =  d~£ 

The  integral  possesses  the  following  properties  that  should  be  taken  into  account 


ou 

J/(#M£¥#  =  o,  iK£)  =  {i,£#2} 


I[fM]  =  0,  fM=n{- 


m 


-)2  exp(- 


y 


) 


(3) 

(4) 

(5) 


2 nkT  2 kT 

The  condition  (4)  leads  to  correct  hydrodynamic  equations,  and  the  fulfillment  of  (5)  by  a  numerical  method 
increases  the  accuracy  of  the  method  in  near  equilibrium  regimes  where  the  solution  is  close  to  a  Maxwellian. 

<>* 

The  collision  integral  at  a  point  q  can  be  presented  in  a  form 


oo  co  2k  bm 


I(£)=\  j  j  f*'~ ff*)gbdbd(pd%d% 


—oo  —oo  0  0 


Let  us  take  £  =  %y  ,  denote  (j)(Zy)  =  S(%  -  %r)  +  S(%,  -  %y)  -  S(%'~  %r)  -  S(%\-  %y)  and  write  the 
integral  in  a  symmetric  fonn 


1 


oo  oo  2k  bm 


Ir  =  -  f  j  j  }  <KZr){f' f.'~  ff.)gbdbdq>d&Z 

—oo  -oo  0  0 


(6) 


For  evaluation  of  the  integral  (6)  consider  a  domain  x  Q  x  2.71  x  bm  in  which  we  introduce  a  uniform 
integration  grid  E,a  dp  >  bv  ■>  <PV  with  Vr  nodes,  such  that  E:u  and  E,^  coincide  with  the  velocity  grid,  and  exclude 
all  values  of  variables  bv,(pv  for  which  post  collision  velocities  ^ ' a  dp  '  fall  outside  ofO  .  Because  the 
points  da  ■>%[]  '  >  in  general  case  don’t  coincide  with  the  velocity  grid  nodes,  a  regularization  of  (6)  is  needed. 

Let  d  an£l  d  denote  the  nearest  vertices  of  cells  inside  which  E, ' a  and  ^  '  lie,  and  d  +s  and  <^/(  _s  denote 

a  pair  of  other  vertices  that  are  located  symmetrically  at  the  opposite  side  of  a  collision  sphere  presented  in  Figure  1 . 
Two  ways  of  the  regularization  define  two  versions  of  the  method  of  evaluation  of  the  collision  integral. 

First  Version:  Conservative  Method  With  Exact  Asymptotic 

Let  us  replace  the  two  last  5  -functions  in  <j){E,y  )  hy  the  expressions 

-£,)  =  (!-  rv )S(dK  -  d )  +  rvS(dv+s  ~  Zr ) . S(£ - d )  =  (1 -  rv )8(dv  ~ )  +  rv5(dv-s  - d ) 
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This  means  that  contributions  to  the  points a  f p  '  are  distributed  between  the  neighboring  nodes.  Denote 
Eo  =  (£  \  )2  +  (£  'A  )2  >■  Ei  =  (£i,  )2  +  (#*,  )2  >  ^2  =  (£v+s)2  +  (^-.s)2  ■  0ne  of  conditions  E]  <  E()  <  E2  ,  or 
E1  <  E0<  Ex  is  true.  The  coefficient  rv  is  defined  by  the  energy  conservation  lawis0  =  (1  —  rv)El  +  rvE., ,  from 
which  it  follows  that  0  <  rv  <  1 .  The  value  of  fa  fp  can  be  defined  by  interpolation.  It  follows  from  the  definition 
of  the  coefficient  rv  that  the  following  formula  is  exact,  when  fp  =  f M  (fp)  at  all  the  grid  nodes 

(7) 

The  integral  (6)  is  evaluated  as  a  sum  simultaneously  in  all  the  nodes  £  of  the  velocity  grid.  With 
B  =  Vxb2m  /(4NC),NC  =  Nv  /  N0,  A„  =  (faJPv  -  faJ,K)gvbv ,  and  Kronnecker’s  symbol  8yp  ,  one  obtains 

Ir  =  +^A)  +  (l-rv)(^  +Sr,My)  +  rv(Sr^+s  +  Sr^)\Av,  (8) 

>-= i 

In  this  method  the  conditions  (4)  and  (5)  are  fulfdled  exactly  for  each  contribution  in  the  sum. 

Second  Version:  Conservative  Method  Without  Interpolation 

Divide  the  integral  (6)  into  two  parts:  the  one  containing  ff„  ,  which  gives  the  contribution  of  “direct”  collisions, 
and  that  with  f  /, ,  which  describes  “inverse”  collisions.  In  the  first  part,  we  apply  the  same  projection  function 
<j)(ff)  .  In  the  second  part,  we  replace  the  kernel  by 

mta,  -£)+(£*  -£)][(i -r;)fKfK+r;fK+j^s}+ 

-£)+(£*  -zrw-r:)fKf»  H8(^+s-^)+8(^v_s-^)y;fK+j^s}gvbv- 

The  coefficient  rv  is  found  from  the  energy  conservation  law  that  gives  rv  =rv/ff  l\rv/Sf]  +(\  —  rv)/s!'f\, 
where  A*!’ =  fp  f ft  ,  and  iff’  =  fA+sfM_s-  Denote  Aj0)  =  fafp  .  The  collision  integral  has  the  form 

f  =  !(<>;,  +  S^m  -  K) A'1'  +  <8?  -  At0)]  +  (8W  +  SM)[(  1  -  rv) A<0)  -  (1  -  rv*)A«] 

v=l 

+(8p,+s,r  +  s^)[rv K0)  -  r;^]}bvgv  (9) 

In  this  version,  the  condition  (5)  is  fulfilled  with  the  accuracy  of  the  order  of  0(h)  for  each  contribution  to  the 
integral  (9).  The  accuracy  for  the  whole  integral  tends  to  0(h2)  when  N v  — »  oo  .  Compared  with  the  first  version, 
the  integral  doesn’t  contain  numerical  diffusion  induced  by  the  interpolation.  Numerical  experiments  have  shown 
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that  for  most  cases,  except  of  the  near  equilibrium  flows,  the  last  method  is  more  efficient.  In  particular,  it  permits 
calculations  with  a  coarse  mesh  in  velocity  space.  The  switching  between  the  two  versions,  if  needed,  is  very  easy. 

Both  methods  are  very  economic  in  terms  of  arithmetic  operations,  because  calculation  of  fv,gv,bv  and  of 
numbers  of  the  functions  fa  ,fp  /„  +s,  f  _s  in  the  array  of  the  distribution  function  can  be  made  in  advance 

and  then  used  for  all  nodes  by  physical  space.  Classical  molecular  potentials  can  be  applied  without  substantial 
slowing  down  of  the  computation,  and  there  is  no  need  for  artificial  molecular  models.  The  8-dimensional 
integration  grid  is  build  by  the  method  [16]  with  an  error  estimate  about  0(N\ versus  0(N^°5)  for  a  Monte 
Carlo  grid. 

Remark  There  is  another  way  to  satisfy  the  conditions  (4)  and  (5):  select  ^ A  and  as  the  pair  of  the  nearest  nodes  on  the 
sphere.  However,  the  approximation  of  real  collisions  in  this  case  is  quite  suspicious. 

Extension  For  Gas  Mixtures  And  Inelastic  Collisions 

A  transformation  of  the  collision  integral  in  impulse  space  (£,£,)  — >  ( p,p* )  is  sufficient  for  extension  of  the 
method  for  a  mixture  of  gases  [10,  11],  Although  the  scheme  of  Figure  1  is  no  more  valid,  the  post  collision 
impulses  pa  ,Pp  lie  symmetrically  against  corresponding  nearest  grid  nodes  pA  ,pfl  ,  and  P^+S,PM  _s,  because 
the  increment  of  the  impulse  can  be  presented  as  8pap  =  (k  +  8 ')/? ,  with  an  integer  vector  k  and  |  8 '  |<  1 . 

* 

The  formulae  for  the  collision  integral  (8)  and  (9),  as  well  as  the  expressions  of  the  splitting  coefficients  r  and  r 
through  the  energies  remain  the  same.  The  realization  of  the  method  for  a  binary  gas  mixture  has  been  done  in  [14]. 
The  case  of  inelastic  collisions  in  a  gas  with  internal  degrees  of  freedom  was  considered  in  [15]  in  the  framework  of 
Wang  Chang  -  Uhlenbeck  kinetic  equation  [17] 

co  2  n  bm 

dfldt  +  Z&ldx^  J  J  { (./;./;  ffJ)lfg,bdbdcpdgj  (10) 

oo  0  0 

Here  Pyl  is  a  probability  of  transition  from  the  energy  levels  (i,  j )  to  ( k ,  l),  and  f.  is  the  velocity  distribution 
function  for  the  level  i.  In  impulse  space,  kinetic  equations  for  binary  reacting  mixtures  have  the  same  form. 

Finite-Difference  Scheme 

The  set  of  N0  equations  (2)  is  solved  by  the  splitting  method  with  a  time  step  T  =  tJ+l  —  t ' 

dfyidt+zrdf;idx= o,  /;■'=//  (id 

dfr/dt  =  ir,  //  f/hl  (i2) 

The  equation  (11)  is  approximated  by  a  second  order  explicit  flux  conservative  scheme.  To  make  the  splitting 
scheme  symmetric,  the  solution  of  (11)  is  repeated  twice  at  a  half- interval  ofr  that  improves  the  accuracy  and 
permits  to  attenuate  the  Courant  condition.  For  the  equation  (12),  different  methods,  including  second  order 
predictor-corrector  scheme,  and  solution  of  an  integral  equation  were  tested.  The  last  scheme  is  written  in  the  form 

t> +1 

f;+l=f;+\mdt 

tj 

By  introducing  the  variable  tv  =  TV  /  Nv  and  an  intermediate  solution  fj~' ,Nv ,  one  obtains  the  scheme 

rj+v!Nv  r  j+(v-l)/Nv  A  y+(v'— 1) / A^v 

J y  ~  J y  +  T  ■  , 

where  A^+Jl  is  the  V  -th  term  of  the  sums  (8)  or  (9).  The  calculation  of  the  collision  integral  is  replaced  by 
continuous  updating  of  the  distribution  function  at  the  relaxation  stage  that  is  convenient  in  realization. 

To  resolve  fast  kinetic  processes,  the  conditions  <  T0,  where  T0  is  the  time  between  collisions,  should  be 
satisfied.  The  size  of  the  mesh  in  configuration  space  can  vary  and  exceed  the  local  mean  free  path  (m.f.p.). 
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EXAMPLES  OF  COMPUTATIONS 


Selected  examples  bellow  illustrate  the  performance  of  the  approach  for  different  ranges  of  physical 
parameters,  rather  than  study  in  details  particular  gas  flows.  To  avoid  unnecessary  complications,  simple  geometries 
of  the  flow  are  considered.  In  all  cases,  boundary  conditions  at  the  surface  are  formulated  as  perfect  accommodation 
with  reflected  Maxwellian  function.  The  computations  are  made  on  a  2.6  Ghz  Pentium-4  PC. 

Near  Equilibrium  and  Low  Speed  Flows 

This  case  is  the  mast  favorable  for  application  of  the  developed  method  by  two  reasons.  First,  the  solution  in  this 
case  is  close  to  a  Maxwellian  function,  and  therefore  its  principal  part  is  evaluated  exactly  in  the  collision  integral. 
Second,  the  needed  area  of  the  velocity  space  is  relatively  small  that  permits  using  a  small  number  of  velocity  grid 
nodes  N0  . 


FIGURE  2.  Low  Speed  Flow  Along  a  Plate 

In  Figure  2,  the  transversal  velocity  in  units  of  1 03  V  is  shown  for  a  hard-sphere  gas  flow  along  a  cold  plate  with 
incoming  Mach  number  M=0.001  and  Knudsen  number  Kn=0.1.  The  plate  is  located  between  X  —  0  andx  =  10, 
in  dimensionless  units  expressed  in  tenns  of  the  m.f.p.  The  velocity  is  in  units  of  ( kT0  /  m)X  2  .  The  velocities  as 

low  as  1 0  5  are  captured  without  problems.  The  computation  until  t  =  1  OOz"0  took  about  15  min  of  CPU  time  and 
usel3  MB  of  computer  memory.  For  this  case,  the  first  method  of  calculation  of  the  collision  integral  had  advantage 
over  the  second  one,  but  for  M>0.02  both  versions  became  equally  efficient. 


FIGURE  3.  Unsteady  Flow  through  a  Quadratic  Orifice 
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Figure  3  presents  a  3D  simulation  of  unsteady  flow  through  a  quadratic  orifice  for  a  density  ratio 
n()  /  n]  =0.995,  where  n J  is  the  density  in  the  left  vessel,  and»0is  that  at  the  right  side  of  the  wall.  The  gas 
temperatures  are  equal  at  both  sides.  The  Knudsen  number  defined  by  the  orifice  side  length  is  Kn=0.2.  Deviations 
of  density  multiplied  by  200  in  the  longitudinal  plane  Z  =  0  are  shown  at  two  moments  of  time  measured  in  units 
of  Z"0  .  The  orifice  is  located  atx  =  0  .  The  problem  is  solved  for  a  gas  with  a  -12  exponential  potential.  The  number 
of  nodes  in  configuration  space  is  12000,  the  number  of  velocity  nodes  is  4224,  computer  memory  usage  is  200 
MB,  and  the  CPU  time  per  r0  is  8  min. 


Hypersonic  Flows 


Consider  now  the  opposite  case  of  big  gradients  and  high  deviations  from  the  thermodynamic  equilibrium  state. 
In  Figure  4,  the  density  and  temperature  fields  for  a  2D  flow  of  a  hard-sphere  gas  are  presented  along  a  cold  plane 
plate  for  M=10  and  Kn=0.01.  The  plate  is  located  along  the  X  axis  with  the  leading  edge  at  .v=0.  All  distances  are 
measured  in  the  units  of  the  m.f.p.  of  the  incoming  flow,  and  flow  parameters  are  related  to  their  unperturbed  values. 


FIGURE  4.  Hypersonic  Flow  along  a  Plate.  The  Density  Field  (Top)  and  the  Temperature  Field  (Bottom) 
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The  computation  is  made  with  5400  nodes  in  configuration  space,  2808  nodes  in  velocity  space,  with  total  memory 
usage  of  about  120  MB,  and  CPU  time  about  3h.  The  test  computation  with  5768  velocity  nodes  gave  practically 
the  same  result  with  maximal  deviations  of  local  heat  and  drag  at  the  plate  of  6%.  In  Figure  5,  the  nonnalized  by 
Wg  /  2  energy  flux  to  the  plate  is  presented  for  two  molecular  models.  For  correct  computation  of  aerodynamic 

reactions,  variable  mesh  in  x  and  y  directions  were  used.  The  mesh  size  along  the  y  axis  near  the  plate  and  along  the 
x  axis  near  the  edges  of  the  plate  were  equal  to  about  0.2  of  the  m.f.p. 

In  Figure  6,  the  density  along  the  symmetry  line  is  shown  for  a  2D  planar  flow  around  an  orthogonal  thin  plate 
with  the  frontal  side  temperature  T f  =  2  and  the  rear  side  temperature  Tr  =  1  for  M=8  and  Kn=0.05.  The  flow  is 

characterized  by  the  fonnation  of  a  sharp  Knudsen  layer  in  front  of  the  plate  and  considerable  rarefaction  behind  the 
plate  and  require  the  application  of  variable  mesh  condensed  near  the  plate.  The  mesh  size  along  the  x  axis  near  the 
plate  was  equal  to  0.05  of  the  m.f.p.  The  stabilization  of  the  solution  behind  the  plate  required  quite  a  long  time  of 
calculation. 

This  example  demonstrates  the  ability  of  the  method  for  computing  flows  with  strong  spatial  variations  of 
hydrodynamic  parameters. 


FIGURE  6.  Hypersonic  Flow  around  an  Orthogonal  Plate.  Density  at  Both  Sides  of  the  Plate  at  .v=0 


Gas  Flows  With  Internal  Degrees  of  Freedom 

As  an  example  of  gas  flows  with  internal  degrees  of  freedom  consider  a  shock  wave  structure  in  Nitrogen  at  room 
temperature  for  M=2,  for  which  it  exists  experimental  data  [18].  The  data  of  rotational  spectrum  and  the  Lennard- 
Jones  potential  for  Nitrogen  were  taken  from  [19,  20],  and  the  R-T  cross  sections  are  from  [21,  22].  In  Figure  7,  the 
normalized  translational,  rotational,  and  longitudinal  temperatures  are  shown.  In  Figure  8,  the  degenerated  rotational 
spectrum  (populations  of  the  levels)  at  5  positions  is  reported,  including  equilibrium  spectrum  before  the  shock  wave 
at  x=-12,  and  behind  the  wave  at  x=l  2.  Figure  9  presents  the  comparison  of  computed  and  experimental  densities. 
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FIGURE  7.  Temperatures 


FIGURE  8.  Rotational  Spectrum 
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FIGURE  9.  Comparison  of  Computed  and  Measured  Densities 
The  computations  were  made  for  35  degenerated  energy  levels  and  require  about  32  MB  of  memory. 
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Appendix:  Output  Of  The  Distribution  Function 


As  a  supplement  for  the  analysis  that  was  made  at  the  level  of  hydrodynamic  parameters,  we  present  the 
distribution  function  for  the  considered  hypersonic  flow  about  the  orthogonal  plate  (contours  at  the  top,  and  the 
surface  at  the  bottom).  Note  the  part  of  the  distribution  that  presents  molecules  going  to  the  plate  (at  the  left),  and  a 
spread  of  the  function  in  t.  direction. 


FIGURE  10.  Distribution  Function  f  (<%x,£y,<!jz,X,y)  atx  =  4 X,y  =  0.25/1 ,  =  h/2 
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